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High energy cosmic ray data collected by the SUGAR Collaboration in the Pilliga State Forest 
have been re-analyzed using up-to-date shower simulations. Complete, inclined angle, Monte Carlo 
simulations reveal 2 events with energies in excess of 7 x 10 19 eV at 95%CL, independent of the 
choice of hadronic interaction model and of chemical composition of the primary. An additional 3 
events, with mean energy > 7 x 10 19 eV, were also re-analyzed in the same manner. A lower bound 
on the flux at the high end of the spectrum, as observed in the Southern sky, is presented on the 
basis of our analysis. 

NUB-3241-TH-03 



INTRODUCTION 

In the mid-60's Greisen, Zatsepin, and Kuzmin (GZK) 
pointed out that if cosmic rays originate at cosmologi- 
cal distances, standard physics implies there would be 
an ultra violet cutoff in the observed spectrum Q . Due 
to a lack of knowledge of the primary species, nowadays 
there is some ambiguity in the definition of the GZK 
energy limit. Gamma-rays produce electromagnetic cas- 
cades and for extragalactic magnetic fields O(nG) the 
high energy component of the spectrum is dramatically 
depleted. The net effect of proton interactions would be 
a pile-up around 5 x 10 19 eV with the spectrum dropping 
sharply thereafter. Nucleus photodisintegration is partic- 
ularly important at energies which excite the giant dipole 
resonance, yielding a cutoff which is species-dependent 
and is slightly shifted to higher energies Q. With this 
in mind, in this work we consider an event to supersede 
the cutoff if its lower energy limit at the 95% CL exceeds 
7 x 10 19 eV. This conforms closely to the most conserva- 
tive criteria delineated in Ref. 0). 

A first hint of a puzzle surfaced in the highest en- 
ergy Fly's Eye event 0] which has no apparent pro- 
genitor within the local supercluster Subsequent 
observations with the Akeno Giant Air Shower Array 
(AGASA) carried strong indication that the cutoff 
was somehow circumvented in the absence of plausible 
nearby sources. Far from adding confirmation to what 
seemed a fascinating discovery, the recent analysis of 
data recorded in monocular mode by the HiRes experi- 
ment has re-infused the field with uncertainty. 

An analysis of the combined data reported by the 
HiRes, the Fly's Eye, and the Yakutsk collaborations is 
supportive of the existence of the GZK cutoff at the > 5a 
(> 3.7a) level. The deviation from GZK depends on the 
set of data used as a basis for power law extrapolation 
from lower energies. In additional input to this analysis, 
it has been recently noted [tj that there may be techni- 
cal problems with the Yakutsk data collection. We have 
nothing further to add to this discussion, but it seems 
worthwhile to analyze any other independent set of data. 



SUGAR DATA REVISITED 



The Sydney University Giant Air-shower Recorder 
(SUGAR) is the largest array ever built in the South- 
ern hemisphere. It ceased operation in February 1979 
after a life of 11 years [[[j- The experiment comprised 
47 detector stations located in the Pilliga State Forest 
(New South Wales, Australia) at latitude 30° 31' S, lon- 
gitude 149° 38' E, and altitude 250 m above sea level. 
Each of these stations consisted of two 6 m 2 conical liq- 
uid scintillator tanks separated by 50 m and buried under 
1.5 ± 0.3 m of earth to shield against the electromag- 
netic component of the shower |ll| . A minimum ionizing 
muon traversing the scintillator produced an average of 
10 photoelectrons, which were detected by photomulti- 
plier tubes (PMTs) looking downward at the scintilla- 
tor. Each tank was equipped with discriminators whose 
thresholds were set to the signal size expected for 3 min- 
imum ionizing muons. A local trigger was generated for 
the case of hits on both tanks within 350 ns of one an- 
other. In this work we are interest in the largest showers, 
which determine the high end of the spectrum (energy 
> 7 x 10 19 eV). In order to obtain a clean sample we 
only analyze events with the shower zenith angle < 45° 
and muon sizes > 10 s particles. Up to 1000 m from the 
core all stations in these showers gave non-zero readings 
and thus well determined densities can be given. We re- 
quire a global trigger of 6 or more stations which did 
record local events, with a minimum of 8 stations trig- 
gered within a time window of 80.5 /is. There are 2 events 
(serial #12420 and #6179) which survive all these cuts. 

A first estimate of the shower arrival direction is made 
by fitting a planar shower front to the station hit times. 
A core search is then carried out using an iterative pro- 
cedure. Next, the lateral distribution function is fitted to 
the observed signals in the stations. In fitting the muon 
density account is taken of both stations which did record 
local events, and those which did not. The SUGAR anal- 
ysis fitted using a Fisher |12| function, which is a modi- 
fied version of the Greisen [HjJ muon structure function 
with the exponent of the distant term being zenith angle 
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Here, N is a normalization constant, 9 is the incident 
zenith angle, r is the perpendicular distance from the 
shower axis, ro = 320 m, a = 0.75, b = 1.50 + 1.86 cos#, 
and 
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In our analysis, we extract the lateral distribution di- 
rectly from Monte Carlo simulations of showers using 
AIRES 2.6.0 0|- For all shower simulations, we spec- 
ify the geographical coordinates of SUGAR to turn on 
the geomagnetic field library in aires [J^. Moreover, 
we discard muons with energies below the muon de- 
tection threshold of a SUGAR scintillator, or (0.75 ± 
0.15) sec 9 GeV. 

A first source of systematic uncertainties is encoded in 
the hadronic interaction model used to process the high 
energy collisions. In the eikonal model 16] of high energy 
hadron-hadron scattering, the unitarized cross section is 
written as 

<?inel = J d 2 b (l-exp{ -2x. oft (s,6)-2x hard (s, £)J) , 

(3) 

where the scattering is compounded as a sum of ladders 
via hard and soft processes through the cikonals x hard and 
X aoft . The leading contenders to approximate the (un- 
known) cross sections at cosmic ray energies, SIBYLL 01 
and QGSjet 0], share the eikonal approximation but 
differ in their ansdtse for the eikonals. In both cases, 
the core of dominant scattering at very high energies is 
the parton-parton minijet cross section. In the QGSjet 
model, this core of the hard eikonal is dressed with a 
soft-pomeron preevolution factor; this amounts to take a 
parton distribution which is Gaussian in the transverse 
coordinate distance \b\. In SIBYLL, the transverse density 
distribution is taken as the Fourier transform of the elec- 
tric form factor, resulting in an exponential (rather than 
Gaussian) falloff of the parton density profile with \b\. 
Some of the shower characteristics resulting from these 
choices are readily predictable: the harder form of the 
SIBYLL form factor allows a greater retention of energy by 
the leading particle, and hence less available for the en- 
suing shower [lflj . This in turn explains the reduction of 
mean secondary multiplicity observed [20j in going from 
QGSjet to SIBYLL simulations. The growth in inelastic 
cross section saturates the log 2 s Froissart bound in both 
cases, but the multiplicative constant turns out to be 
larger for SIBYLL, enhancing the cross section by about 
30% at energies of 10 19-20 eV. However, the shower char- 
acteristics are the result of exponentially many hadronic 



collisions in air, and for these the larger SIBYLL cross sec- 
tion is somewhat offset by its smaller multiplicity, reduc- 
ing the differences between the two schemes down to the 
level of 10-20%. These and other properties of the two 
schemes have been amply displayed in a recent work |2l| . 

Eight stations were involved in the event #12420, all 
with local triggers (i.e., 16 scintillators). A map of the 
recorded muon densities in a plane perpendicular to the 
shower axis was reported by the SUGAR Collaboration in 
Ref. jl^l • The north tank of one of the stations recorded 
a saturated response (density > 800 m~ 2 ). The standard 
analysis of the SUGAR Collaboration placed the shower 
core on this tank. However, the collaboration concluded 
that in such a big shower saturated densities may occur 
at distances up to several hundred meters from the core. 
Moreover, for large signals, afterpulsing in the PMTs can 
lead to an overestimate of the muon density. The loca- 
tion of the core was best re-estimated by assigning lower 
energy densities to this station. Since the probability 
of recording a saturated response from an actual den- 
sity less than 100 m~ 2 is negligible, the SUGAR Col- 
laboration concluded that this value of density in both 
tanks would yield a minimum value on the primary en- 
ergy. In our analysis we adopt this conservative lower 
bound. Note that this second source of systematic error 
can lead to an underestimate of the energy. The best x 2 
fit for QGSjetOl (sibyll 2.1), assuming a proton primary, 
is given in Fig. ^ (Fig. El- Note that aires + sibyll sim- 
ulations result in an energy which exceeded the AIRES + 
QGSjet value by about 10%. This agrees with the analysis 
of the AGASA Collaboration [H . 

The energy that we find for event#12420 using QGSjet 
is about 5% above the energy reported by the SUGAR 
Collaboration llfl. their estimate being based on the 
Hillas E model 123 
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In further exploring the energy systematics, we have gen- 
erated proton vertical showers using AiRES+QGSjet and 
compared (in Fig. |3J) the A^ crt — E relation with that 
given in Eq. (0J. As can be seen from Fig. El this conver- 
sion formula underestimates the energy by about 25%, 
for a given -/V™ rt , when compared with our simulations. 
However, as noted above, the overall underestimate of 
the SUGAR assignment for the energy from our best fit 
is about 5%. This is readily explained by noting that the 
normalization constant N (used by the SUGAR Collab- 
oration as the muon size) and obtained through the fit to 
the Fisher function overestimates the number of muons 
(as found through our simulation) by about 17% (see Ta- 
bleP). Since the equivalent vertical muon size and energy 
bear a near-linear relationship to one another (see Fig.|2Jl, 
these systematic effects nearly offset one another, bring- 
ing the SUGAR energy estimate near to that obtained 
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FIG. 1: Best fit to the muon density (*) for event #12420 
using AIRES + QGSjet. The estimated primary energy is 
1.1 x 10 20 eV, with xVd.o.f. = 9.4/7. The recorded par- 
ticle densities as reported by the SUGAR Collaboration are 
indicated with squares. The solid line indicates the best fit 
reported by SUGAR Collaboration using the Fisher function 
( x 2 /d.o.f. = 10.3/7). 



FIG. 3: Total number of muons as a function of the en- 
ergy shower. The simulations were carried out with AIRES 
+ QGSjet, assuming protons as primaries. The particles were 
injected vertically at the top of the atmosphere and the detec- 
tor was placed at 250 m above sea level. Muons with energies 
below the SUGAR threshold (0.75 GeV) are not taken into 
account in the simulations. The dashed line indicates the 
Hillas parametrization (see main text). 




FIG. 2: Best fit to the muon density (*) for event #12420 us- 
ing AIRES + SIBYLL. The estimated energy is 1.2 x 10 20 eV, 
with xVd.o.f. = 9.1/7. The recorded particle densities 
as reported by the SUGAR Collaboration are indicated 
with squares. The solid line indicates the best fit re- 
ported by SUGAR Collaboration using the Fisher function 
(x 2 /d.o.f. = 10.3/7). 



with AiRES+QGSjet. Implicitly, this discussion supports 
the consistency of the constant intensity cut method of 
mapping muon sizes from inclined to vertical |25|. 

A third source of systematic uncertainty is hidden in 
the chemical composition. As a first approximation, a 
shower produced by a nucleus with energy E A and mass 
A can be modeled by the collection of A proton showers, 
each with A^ 1 of the nucleus energy. Now, using Eq. @ a 
straightforward calculation shows that the total number 
of muons produced by the superposition of A individual 
proton showers is, NY\ CIt cx A(E A /A) - 93 . Consequently, 
one expect a cosmic ray nucleus to produce about A 0,07 
more muons than a proton. Conversely, a given muon 
size will imply a smaller energy for a complex nucleus. 
Thus, we can estimate the lower energy limit of the event 
#12420, by simulating several sets of showers for 56 Fe 
primaries [26j. The QGSjet package was used to process 
the high energy hadronic interactions. Energies below 
7.0 x 10 19 eV are excluded at the 95% CL. 

The analysis of the second event (#6179) is less 
straightforward. This is because there were stations in 
which one of the tanks did not fire [j^. In assigning 
the trigger intensity of these tanks, the SUGAR Col- 
laboration took as a maximum the density in the firing 
tanks, and as a minimum the density averaged over all 
the tanks that were expected to fire. Data scaled to a 
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TABLE I: Primary energy in units of 10 eV for proton pri- 
maries estimated using AIRES + QGSjet. The corresponding 
muon sizes and normalization constant of the Fisher function 
are also given in units of 10 8 particles. The last 3 columns 
indicate the incident zenith angle and the arrival directions 
in equatorial coordinates (a, 8) with respect to the reference 
system B1950 fiH . The estimated angular uncertainty for all 
these events is 4.3°. 



Number 


M 


E 




[deg.] 


a [deg.] 


S [deg.] 


12420 


3.47 


1.10 


2.99 ±0.35 


43 


147.0 


-43.6 


6179 


4.06 


1.03 


3.56 ±0.45 


24 


121.5 


-32.0 


7329 


2.63 


0.77 


2.24 ±0.34 


38 


135.0 


04.3 


1696 


2.55 


0.71 


2.20 ±0.27 


33 


331.5 


-32.2 


6239 


2.27 


0.70 


1.96 ±0.27 


40 


280.5 


-55.6 



muon size of 4.5 x 10 s is shown in Ref. 22] . Because 
of the resulting decrease in the quality of the statistics, 
the likelihood in this case did not converge reliably to 
a global minimum. However, a 95% confidence interval, 
7.0 x 10 19 eV < E < 1.2 X 10 20 eV, was readily obtained. 

For a more extended comparison with the analysis re- 
ported by the SUGAR Collaboration, we relax the global 
trigger designated above, and follow several more recent 
analyses of SUGAR data |22j by considering events in 
which 5 pairs of scintillators fired within a time window of 
80.5 /xs. For this kind of global trigger, the uncertainty 
in the arrival direction is manageable at an average of 
4.3° [Hf. (For 3- and 4- fold triggers, this rises beyond 10 
degrees.) With this new global trigger, we simulate show- 
ers with AIRES + QGSjet and find that by increasing the 
energy reported by the SUGAR Collaboration by about 
5% we reasonably conform to the lateral profiles of the 
events as described by the Fisher function, with normal- 
ization constants taken from |lfj| and given in Table [I] 
As a result, we find 4 events (including #6179 discussed 
above) with energies > 7 x 10 19 eV. These qualitative fits 
are shown in Fig. and the parameters characterizing 
these events are detailed in Table H] 

A closer examination of Fig. 0] reveals a systematic de- 
viation between the QGSjet and Fisher function lateral 
distribution profiles. As can be seen in Fig. [SJ the di- 
vergences become more severe when SIBYLL is used to 
process the hadronic interactions. Of course, this com- 
parison with the Fisher function plays no role in discrim- 
inating between these schemes. 

ENERGY SPECTRUM 

The effective detecting area of the array at any time 
is given by the product of the active geometric surface 
area of the array, S, times the probability p(Af, 0) that a 
shower falling within S will be detected, 




FIG. 4: Lateral distributions obtained using air-shower sim- 
ulations with AIRES ± QGSjet (*) superimposed over the best 
fits of the Fisher function (solid line) as reported by the 
SUGAR Collaboration. The corresponding primary energies 
used in the simulations are those given in Table 




FIG. 5: Lateral distributions obtained using air-shower simu- 
lations with AIRES ± SIBYLL (*) superimposed over the best 
fits of the Fisher function (solid line) as reported by the 
SUGAR Collaboration. The corresponding primary energies 
in units of 10 19 eV used in the simulations are: 11.3, 8.5, 7.8, 
and 7.7, for #6179, #7329, #1696 and #6239, respectively. 



A eS (M,e,t) = S(t)p(M,6). (5) 
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The SUGAR Collaboration reported [23 an angle specific 
exposure for showers with AT = 4 x 10 s and 9 = 32°, 



26 



£(32 c 



dtS{t)p{4 x 10 8 ,32°) = 540 km 2 yr, (6) 



where T = 11 yr. This value of p was obtained by allow- 
ing a global trigger of 3 or more stations. Because the 
detection probability decreases with smaller muon size 
(as is the case for the events in Table and with our re- 
striction to 5 triggers, Eq. © leads to an upper limit on 
the integrated exposure for the parameter space relevant 
to this study, 



£(9 < 45 c 



2tt £(32°) 



V2/2 P(32°) 



cos 9 dcos9 



862 km sr yr, 



(7) 



where p(9) oc cos (9 is the probability of detection at 
zenith angle 9 pffll ]. 

The observed number of events in an energy bin A is 
given by 



N 



2.3 £ EJ(E) Alog 10 £; 



so that 



E 3 J(E) 



E 2 



A 



2.3 £ A\og 10 E 



(8) 



(9) 



where E is the value of E at the center of the logarithmic 
interval. 

Following common practice, we adopt a bin size 
Alog 10 i? = 0.1. As previously discussed, we have 
adopted a lower energy cut of 10 19 - 85 eV. The events 
listed in Tableware then grouped in two bins, and the re- 
sulting contributions to the cosmic ray spectrum as given 
by Eq. © are shown in Fig. [B| The error on the mean 



Alog 10 [£ J(E)} = Alog 10 iV 



0.343 A In AT 
AA 

= 0.343 — 

A 



is taken as Poisson (i.e., A A = v A). 

On the same figure are shown the data points for the 
flux as reported by the AGASA 23], Fly's Eye H3, and 
HiRes 7] collaborations. Also shown is the flux obtained 
by a recalibration of the 4 highest energy events observed 
at the Haverah Park facility [23] • We note that the 
AGASA flux plotted is taken from the most recent analy- 
sis published by the AGASA Collaboration [2^] , contain- 
ing events with zenith angle < 45°, recorded until July 
2002. The energy is calibrated using an average over 
hadronic interaction models and primary species (p and 
56 Fe). It is of interest to the present work that the ener- 
gies obtained in this manner differ by only about 2% with 
respect to those obtained via aires + QGSjet simulation. 
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FIG. 6: Upper end of the cosmic ray energy spectrum as 
observed by 5 different experiments (for details see main text). 



SUMMARY AND DISCUSSION 

(1) We have re-examined cosmic ray data taken over 
two decades ago at the SUGAR facility, using up-to-date 
simulation methods to describe shower evolution in the 
atmosphere, with resulting correlation between muon size 
and energy. We conclude that there is strong evidence 
for the existence of events superseding the GZK cutoff. 
Specifically, there are two events in the sample with en- 
ergy > 7x 10 19 eV at the 95%CL. This statement includes 
consideration of systematic uncertainties stemming from 
choice of hadronic interaction model and chemical com- 
position. 

(2) As part of our exploration of the systematics, we 
have compared our study with methods used by the 
SUGAR Collaboration for energy assignment. We found 
that through a compensation of errors in assessing muon 
size from fits to the Fisher function and the normaliza- 
tion of the Hillas relation jiij , the final SUGAR energies 
are consistently about 5% below those found via aires 
+ QGSjet simulation. Remarkably, we find (Fig. |3J) that 
the exponent in the power law A™ rt = JCE a is consistent 
with Hillas's a — 0.93, but with a shift in the constant JC 
upward by 25%. A detailed recent study j^lj shows that 
in fact the exponent a is more properly taken to have a 
logarithmic energy-dependence. 

(3) Using a conservative estimate on the acceptance 
of the SUGAR facility, we have been able to present 
lower bounds on the cosmic ray flux in the case of a 
proton primary with QGSjet as the hadronic interaction 
model. This lower bound would prevail for SIBYLL inde- 
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pendent of the primary composition. The lower bounds 
would be slightly modified in the case that QGSjet suc- 
cessfully models the hadronic interaction at high energies 
and the primary composition is dominated by heavy nu- 
clei. Within a standard deviation, our data points are 
consistent with both AG AS A and HiRes. Because of the 
Southern hemisphere exposure of SUGAR, this comprises 
interesting support for North-South isotropy in the ar- 
rival direction of the highest energy cosmic rays. 

(4) An extrapolation of our analysis to the entire 
SUGAR data set would involve observations using a 3 
station global trigger, with its attendant loss of angular 
resolution beyond 10°. Because of the similar procedures 
for energy assignment in ground arrays, the consistency 
of the SUGAR spectrum with AGASA observations sup- 
ports recent anisotropy studies 33] of low-multipole in- 
homogeneities that combine data taken from these 2 ex- 
periments. 

(5) To give a visual impression of the Southern cos- 
mic ray sky, in Fig.Qwe show the arrival direction of the 
events given in Table[I]together with the potential sources 
included in the first budget of the anisotropy search pro- 
gram for the Pierre Auger Observatory |3J| . 

(6) As part of our study, we have noted an interest- 
ing differentiating signature between the QGSJET and 
SIBYLL models. Previous work [2(| has shown that for a 
fixed (proton) primary energy, muon densities obtained 
in shower simulations using SIBYLL fall more rapidly with 
lateral distance to the shower core than those obtained 
using QGSJET. This can be understood as a manifestation 
of the enhanced leading particle effect in SIBYLL, which 
can be traced to the relative hardness of the electromag- 
netic form factor profile function. From Figs. |@} and JSJ, 
it is apparent that even when the energy is left as a free 
parameter, the curvature of the distribution (dr p^/dr 2 ) 
is measurably different in the two cases, and could possi- 
bly serve as a discriminator between hadronic interaction 
models with sufficient statistics. 

In summary, the footprints analyzed in this paper are 
an interesting addition to the data set, but not enough 
to weigh in on one side or the other with respect to the 
GZK question. The Pierre Auger Observatory now 
coming on line with its superior energy resolution and 
statistics will supply the final verdict on the GZK cutoff. 
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FIG. 7: Arrival direction of the events given in Table|H(o) in 
equatorial coordinates B1950. Also shown in the figure are 
the position in the sky of potential cosmic ray sources: Cen A 
(*) H, NGC3256 (■) Hi, and NGC253 (•) || The solid 
line indicates the position of the Galactic plane. 
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